R Notebook
GARCH in R
install.packages("PerformanceAnalytics", type = "source")## Installing package into '/opt/homebrew/lib/R/4.4/site-library'
## (as 'lib' is unspecified)library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(zoo)
library(PerformanceAnalytics)##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendgetwd()## [1] "/Users/michal/Documents/Code/metals"# setwd("C:/Users/user/OneDrive/Dokumenty/R/classes2024")
list.files()## [1] "__pycache__" "bloomberg_data"
## [3] "combined_data.csv" "data"
## [5] "data_loading.ipynb" "data_vis.ipynb"
## [7] "fetch_dailymetalprice.py" "GARCH_S&P500.ipynb"
## [9] "GARCH.html" "GARCH.ipynb"
## [11] "GARCH.Rmd" "README.md"
## [13] "sp500prCl.csv" "utils"setwd("/Users/michal/Documents/Code/metals")
sp500prices <- read.csv("sp500prCl.csv", dec = ".", sep = ";")
# sp500prices
# Plot daily S&P 500 prices
plot(sp500prices$Price)
class(sp500prices)## [1] "data.frame"stockdata <- sp500prices
date <- as.POSIXct(strptime(stockdata$Data, format = "%Y-%m-%d")) # "%Y-%m-%d" "%d.%m.%Y"
# date
sp500prices1 <- xts(stockdata[2:ncol(stockdata)], order.by = date) #
rm(stockdata)
head(sp500prices1$Price)## Price
## 1989-01-03 275.31
## 1989-01-04 279.43
## 1989-01-05 280.01
## 1989-01-06 280.67
## 1989-01-09 280.98
## 1989-01-10 280.38# Compute daily returns
# sp500ret <- diff(log(sp500prices1))
sp500ret <- CalculateReturns(sp500prices1) # watch the dog! first ret is zero!
sp500ret <- sp500ret[-1, ]
# Check the class of sp500ret
class(sp500ret)## [1] "xts" "zoo"# Plot daily returns
plot(sp500ret)
# Compute the daily standard deviation for the complete sample
sd(sp500ret)## [1] 0.0109947# Compute the annualized volatility for the complete sample
sqrt(252) * sd(sp500ret)## [1] 0.1745354# Compute the annualized standard deviation for the year 2009
sqrt(252) * sd(sp500ret["2009"])## [1] 0.2728327# Compute the annualized standard deviation for the year 2017
sqrt(252) * sd(sp500ret["2017"])## [1] 0.06685658# Showing two plots on the same figure
par(mfrow = c(2, 1))
# Compute the rolling 1 month estimate of annualized volatility
chart.RollingPerformance(
R = sp500ret["2000::2017"], width = 22,
FUN = "sd.annualized", scale = 252, main = "One month rolling volatility"
)
# Compute the rolling 3 months estimate of annualized volatility
chart.RollingPerformance(
R = sp500ret["2000::2017"], width = 66,
FUN = "sd.annualized", scale = 252, main = "Three months rolling volatility"
)
## Compute absolute value of the prediction errors
# Compute the mean daily return
m <- mean(sp500ret)
# Define the series of prediction errors
e <- sp500ret - m
# Plot the absolute value of the prediction errors
par(mfrow = c(2, 1), mar = c(3, 2, 2, 2))
plot(abs(e))
# Plot the acf of the absolute prediction errors
acf(abs(e), lag.max = 20)
# GARCH:
# Define parameters and numerics
alpha <- 0.1
beta <- 0.8
omega <- var(sp500ret) * (1 - alpha - beta)
e <- sp500ret - mean(sp500ret)
e2 <- e^2
nobs <- length(sp500ret)
predvar <- rep(NA, nobs)
# Compute the predicted variances
predvar[1] <- var(sp500ret)
for (t in 2:nobs) {
predvar[t] <- omega + alpha * e2[t - 1] + beta * predvar[t - 1]
}
# Create annualized predicted volatility
ann_predvol <- xts(sqrt(252) * sqrt(predvar), order.by = time(sp500ret))
# Plot the annual predicted volatility in 2008 and 2009
plot(ann_predvol["2008::2009"], main = "Ann. S&P 500 vol in 2008-2009")
# Time for rugarch package
# install.packages("rugarch")
library(rugarch)## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:stats':
##
## sigma# Specify a standard GARCH model with constant mean
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "sGARCH"),
distribution.model = "norm"
)
# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
# See the resutls:
garchfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000572 0.000087 6.59621 0.000000
## omega 0.000001 0.000002 0.67008 0.502805
## alpha1 0.078259 0.027893 2.80567 0.005021
## beta1 0.910553 0.029458 30.91016 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000572 0.000853 0.670332 0.50265
## omega 0.000001 0.000058 0.021204 0.98308
## alpha1 0.078259 0.901416 0.086818 0.93082
## beta1 0.910553 0.946490 0.962031 0.33603
##
## LogLikelihood : 24073.8
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.5881
## Bayes -6.5844
## Shibata -6.5881
## Hannan-Quinn -6.5869
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.311 0.25213
## Lag[2*(p+q)+(p+q)-1][2] 1.582 0.34264
## Lag[4*(p+q)+(p+q)-1][5] 7.253 0.04502
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.07822 0.7797
## Lag[2*(p+q)+(p+q)-1][5] 5.15209 0.1414
## Lag[4*(p+q)+(p+q)-1][9] 6.34246 0.2607
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1801 0.500 2.000 0.6713
## ARCH Lag[5] 0.2836 1.440 1.667 0.9446
## ARCH Lag[7] 0.8098 2.315 1.543 0.9424
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 346.8562
## Individual Statistics:
## mu 0.04305
## omega 30.47592
## alpha1 0.20103
## beta1 0.29761
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 2.6152 8.935e-03 ***
## Negative Sign Bias 0.3192 7.496e-01
## Positive Sign Bias 2.3649 1.806e-02 **
## Joint Effect 33.3642 2.698e-07 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 223.6 6.447e-37
## 2 30 232.5 1.196e-33
## 3 40 270.2 2.282e-36
## 4 50 293.3 1.537e-36
##
##
## Elapsed time : 0.2587059# Use the method sigma to retrieve the estimated volatilities
garchvol <- sigma(garchfit)
garch_mean <- fitted(garchfit)
# Plot the volatility for 2017
plot(garchvol["2017"])
# Forecasts:
# Compute unconditional volatility
sqrt(uncvariance(garchfit))## [1] 0.01051693# Print last 10 ones in garchvol
tail(garchvol, 10)## [,1]
## 2017-12-15 0.005071375
## 2017-12-18 0.005493771
## 2017-12-19 0.005524124
## 2017-12-20 0.005491357
## 2017-12-21 0.005371077
## 2017-12-22 0.005259472
## 2017-12-26 0.005148613
## 2017-12-27 0.005057916
## 2017-12-28 0.004953331
## 2017-12-29 0.004868582# Forecast volatility 5 days ahead and add
garchforecast <- ugarchforecast(fitORspec = garchfit, n.ahead = 5)
# Extract the predicted volatilities and print them
print(sigma(garchforecast))## 2017-12-29
## T+1 0.005041038
## T+2 0.005134710
## T+3 0.005225683
## T+4 0.005314106
## T+5 0.005400117Volatility targeting
# When you target a portfolio with 5% annualized volatility, and the annualized volatility of the risky asset is, then you should invest in the risky asset.
# Compute the annualized volatility
annualvol <- sqrt(252) * sigma(garchfit)
# Compute the 5% vol target weights
vt_weights <- 0.05 / annualvol
# Compare the annualized volatility to the portfolio weights in a plot
plot(merge(annualvol, vt_weights), multi.panel = TRUE)
# Application of the GARCH models to the data:
# generate ret series"
# list(mu = 0, ar1 = 0, ma1 = 0, omega = 6*10^(-7), alpha1 = 0.07, beta1 = 0.9, skew = 0.9, shape = 5)
ret <- sp500ret
# Plot the return series
plot(ret)
# Specify the garch model to be used
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "sGARCH"),
distribution.model = "sstd"
)
# Estimate the model
garchfit <- ugarchfit(data = ret, spec = garchspec)
garchvol <- sigma(garchfit)
garch_mean <- fitted(garchfit)# Inspect the coefficients
coef(garchfit)## mu omega alpha1 beta1 skew shape
## 5.648052e-04 6.398380e-07 7.521339e-02 9.218343e-01 9.435684e-01 6.315897e+00# See the whole output
garchfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000565 0.000087 6.4565 0.000000
## omega 0.000001 0.000000 2.7701 0.005604
## alpha1 0.075213 0.004718 15.9414 0.000000
## beta1 0.921834 0.004427 208.2420 0.000000
## skew 0.943568 0.014671 64.3136 0.000000
## shape 6.315897 0.454050 13.9101 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000565 0.000081 6.9867 0.000000
## omega 0.000001 0.000001 1.0855 0.277682
## alpha1 0.075213 0.020742 3.6262 0.000288
## beta1 0.921834 0.018987 48.5506 0.000000
## skew 0.943568 0.014018 67.3091 0.000000
## shape 6.315897 0.627827 10.0599 0.000000
##
## LogLikelihood : 24279.25
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.6438
## Bayes -6.6382
## Shibata -6.6438
## Hannan-Quinn -6.6419
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.197 0.27384
## Lag[2*(p+q)+(p+q)-1][2] 1.498 0.36140
## Lag[4*(p+q)+(p+q)-1][5] 7.382 0.04186
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.002836 0.9575
## Lag[2*(p+q)+(p+q)-1][5] 4.322017 0.2165
## Lag[4*(p+q)+(p+q)-1][9] 5.564033 0.3514
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.2192 0.500 2.000 0.6396
## ARCH Lag[5] 0.2926 1.440 1.667 0.9422
## ARCH Lag[7] 1.0704 2.315 1.543 0.9018
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 748.9782
## Individual Statistics:
## mu 0.07006
## omega 156.38935
## alpha1 0.37783
## beta1 0.56251
## skew 0.41544
## shape 0.71577
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 2.7461 6.045e-03 ***
## Negative Sign Bias 0.1844 8.537e-01
## Positive Sign Bias 2.5478 1.086e-02 **
## Joint Effect 33.9210 2.059e-07 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 53.96 3.349e-05
## 2 30 61.00 4.632e-04
## 3 40 71.14 1.263e-03
## 4 50 83.73 1.463e-03
##
##
## Elapsed time : 0.38465# Compute the standardized returns
stdret <- residuals(garchfit, standardize = TRUE)
# Compute the standardized returns using fitted() and sigma()
stdret <- (ret - fitted(garchfit)) / sigma(garchfit)
chart.Histogram(stdret,
methods = c("add.normal", "add.density"),
colorset = c("gray", "red", "blue")
)
Leverage effect:
# Specify the GJR GARCH model
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "sstd"
)
# Estimate the model and compute volatility
gjrgarchfit <- ugarchfit(data = sp500ret, spec = garchspec)
gjrgarchvol <- sigma(gjrgarchfit)Other specifications:
# Specify AR(1)-GJR GARCH model
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(1, 2)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "sstd"
)
# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
# Print the first four coefficients
coef(garchfit)[c(1:4)]## mu ar1 ma1 ma2
## 0.0004002116 0.6689868348 -0.7027022563 -0.0109852393APARCH
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 1)),
variance.model = list(model = "apARCH"),
distribution.model = "std"
)
# Estimate the model
aparchfit <- ugarchfit(data = sp500ret, spec = garchspec, solver = "hybrid")
aparchvol <- sigma(aparchfit)
aparchmean <- fitted(aparchfit)A remark on the skewed Student distribution:
The conditional density to use for the innovations. Valid choices are “norm” for the normal distibution, “snorm” for the skew-normal distribution, “std” for the student-t, “sstd” for the skew-student, “ged” for the generalized error distribution, “sged” for the skew-generalized error distribution, “nig” for the normal inverse gaussian distribution, “ghyp” for the Generalized Hyperbolic, and “jsu” for Johnson’s SU distribution.
# Compare volatility
plotvol <- plot(abs(sp500ret), col = "grey")
plotvol <- addSeries(gjrgarchvol, col = "red", on = 1)
plotvol <- addSeries(garchvol, col = "blue", on = 1)
plotvol <- addSeries(aparchvol, col = "green3", on = 1)
plotvol
GARCH-in-Mean specification and estimation
gim_garchspec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), archm = TRUE, archpow = 2), variance.model = list(model = "gjrGARCH"), distribution.model = "sstd")
gim_garchfit <- ugarchfit(data = sp500ret, spec = gim_garchspec)
gim_garchfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000217 0.000073 2.9707 0.002972
## archm 2.038474 0.544626 3.7429 0.000182
## omega 0.000001 0.000002 0.5719 0.567390
## alpha1 0.000740 0.004539 0.1631 0.870441
## beta1 0.913129 0.029509 30.9445 0.000000
## gamma1 0.148437 0.058562 2.5347 0.011255
## skew 0.922223 0.012843 71.8092 0.000000
## shape 7.018567 0.902143 7.7799 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000217 0.001314 0.165245 0.86875
## archm 2.038474 23.768012 0.085765 0.93165
## omega 0.000001 0.000045 0.027598 0.97798
## alpha1 0.000740 0.023691 0.031248 0.97507
## beta1 0.913129 0.620573 1.471430 0.14117
## gamma1 0.148437 1.243053 0.119413 0.90495
## skew 0.922223 0.145021 6.359238 0.00000
## shape 7.018567 22.134566 0.317086 0.75118
##
## LogLikelihood : 24374.79
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.6694
## Bayes -6.6619
## Shibata -6.6694
## Hannan-Quinn -6.6668
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.3488 0.5548
## Lag[2*(p+q)+(p+q)-1][2] 0.4322 0.7263
## Lag[4*(p+q)+(p+q)-1][5] 4.6044 0.1877
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 3.022 0.08214
## Lag[2*(p+q)+(p+q)-1][5] 4.002 0.25384
## Lag[4*(p+q)+(p+q)-1][9] 5.162 0.40594
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.359 0.500 2.000 0.2437
## ARCH Lag[5] 1.546 1.440 1.667 0.5803
## ARCH Lag[7] 2.579 2.315 1.543 0.5963
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 397.6714
## Individual Statistics:
## mu 0.2736
## archm 0.1413
## omega 59.5030
## alpha1 1.2735
## beta1 1.3664
## gamma1 0.5872
## skew 0.8253
## shape 0.9278
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.89 2.11 2.59
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 3.234 1.228e-03 ***
## Negative Sign Bias 2.795 5.198e-03 ***
## Positive Sign Bias 2.180 2.929e-02 **
## Joint Effect 28.717 2.568e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 52.19 6.203e-05
## 2 30 54.06 3.184e-03
## 3 40 66.23 4.188e-03
## 4 50 69.03 3.113e-02
##
##
## Elapsed time : 0.882395# Predicted mean returns and volatility of GARCH-in-mean
gim_mean <- fitted(gim_garchfit)
gim_vol <- sigma(gim_garchfit)The comparison of different conditional means and volatilites
# Correlation between predicted return using one of the previous models and GARCH-in-mean models
cor(aparchmean, gim_mean)## [,1]
## [1,] 0.09720179cor(garch_mean, gim_mean)## [,1]
## [1,] -0.1126776# Correlation between predicted volatilities across models
cor(merge(garchvol, gjrgarchvol, aparchvol, gim_vol))## garchvol gjrgarchvol aparchvol gim_vol
## garchvol 1.0000000 0.9617339 0.9773612 0.9623802
## gjrgarchvol 0.9617339 1.0000000 0.9771208 0.9995686
## aparchvol 0.9773612 0.9771208 1.0000000 0.9773793
## gim_vol 0.9623802 0.9995686 0.9773793 1.0000000Avoid unneccessary complexity:
- restrict parameters (setbounds(), setfixed())
# Print the flexible GARCH parameters # flexgarchfit
print(coef(gim_garchfit))## mu archm omega alpha1 beta1 gamma1
## 2.170667e-04 2.038474e+00 1.252247e-06 7.403016e-04 9.131289e-01 1.484370e-01
## skew shape
## 9.222228e-01 7.018567e+00# Restrict the flexible GARCH model by impose a fixed ar1 and skew parameter
rflexgarchspec <- gim_garchspec
setfixed(rflexgarchspec) <- list(archm = 1, skew = 1)
# Estimate the restricted GARCH model #EURUSDret
rflexgarchfit <- ugarchfit(data = sp500ret, spec = rflexgarchspec)
rflexgarchfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000419 0.000081 5.17139 0.00000
## archm 1.000000 NA NA NA
## omega 0.000001 0.000001 1.48500 0.13754
## alpha1 0.001013 0.003950 0.25648 0.79758
## beta1 0.916079 0.010846 84.46270 0.00000
## gamma1 0.144571 0.014695 9.83783 0.00000
## skew 1.000000 NA NA NA
## shape 6.719495 0.341236 19.69165 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000419 0.000111 3.761047 0.000169
## archm 1.000000 NA NA NA
## omega 0.000001 0.000006 0.190751 0.848721
## alpha1 0.001013 0.024118 0.042008 0.966492
## beta1 0.916079 0.093999 9.745584 0.000000
## gamma1 0.144571 0.149394 0.967716 0.333186
## skew 1.000000 NA NA NA
## shape 6.719495 3.063479 2.193419 0.028277
##
## LogLikelihood : 24360.86
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.6662
## Bayes -6.6605
## Shibata -6.6662
## Hannan-Quinn -6.6642
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5871 0.4435
## Lag[2*(p+q)+(p+q)-1][2] 0.7641 0.5815
## Lag[4*(p+q)+(p+q)-1][5] 5.7898 0.1008
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.680 0.1016
## Lag[2*(p+q)+(p+q)-1][5] 3.700 0.2939
## Lag[4*(p+q)+(p+q)-1][9] 4.884 0.4466
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.334 0.500 2.000 0.2480
## ARCH Lag[5] 1.557 1.440 1.667 0.5774
## ARCH Lag[7] 2.608 2.315 1.543 0.5905
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 430.727
## Individual Statistics:
## mu 0.2277
## omega 78.1497
## alpha1 1.1330
## beta1 1.2566
## gamma1 0.6079
## shape 0.9506
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 3.087 2.029e-03 ***
## Negative Sign Bias 2.544 1.099e-02 **
## Positive Sign Bias 2.168 3.019e-02 **
## Joint Effect 26.943 6.050e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 64.33 7.829e-07
## 2 30 73.14 1.113e-05
## 3 40 88.10 1.164e-05
## 4 50 101.40 1.612e-05
##
##
## Elapsed time : 0.7889161# Compare the volatility of the unrestricted and restriced GARCH models
plotvol <- plot(abs(sp500ret), col = "grey")
plotvol <- addSeries(sigma(gim_garchfit), col = "black", lwd = 4, on = 1)
plotvol <- addSeries(sigma(rflexgarchfit), col = "red", on = 1)
plotvol
Another solution: setbounds()
# Define bflexgarchspec as the bound constrained version
bflexgarchspec <- gim_garchspec
setbounds(bflexgarchspec) <- list(alpha1 = c(0.05, 0.2), beta1 = c(0.7, 0.95))
# Estimate the bound constrained model
bflexgarchfit <- ugarchfit(data = sp500ret, spec = bflexgarchspec)
bflexgarchfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : sstd
##
## Convergence Problem:
## Solver Message:# Inspect coefficients
coef(bflexgarchfit)## NULL# Compare forecasts for the next ten days:
forecast_gim_garch <- ugarchforecast(gim_garchfit, n.ahead = 10)
series_forecast_gim <- forecast_gim_garch@forecast$seriesFor
sigma_forecast_gim <- forecast_gim_garch@forecast$sigmaFor
forecast_garch <- ugarchforecast(garchfit, n.ahead = 10)
sigma_forecast_garch <- forecast_garch@forecast$sigmaFor
cbind(sigma_forecast_gim, sigma_forecast_garch)## 2017-12-29 2017-12-29
## T+1 0.004782920 0.004617103
## T+2 0.004878952 0.004697538
## T+3 0.004971807 0.004775718
## T+4 0.005061678 0.004851765
## T+5 0.005148739 0.004925786
## T+6 0.005233148 0.004997880
## T+7 0.005315048 0.005068141
## T+8 0.005394570 0.005136651
## T+9 0.005471833 0.005203489
## T+10 0.005546948 0.005268728#Variance targeting Financial return volatility clusters through time: periods of above average volatility are followed by period of below average volatility. The long run prediction is that:
when volatility is high, it will decrease and revert to its long run average.
when volatility is low, it will increase and revert to its long run average.In the estimation of GARCH models we can exploit this mean reversion behavior of volatility by means of volatility targeting. We then estimate the GARCH parameters in such a way that the long run volatility implied by the GARCH model equals the sample standard deviation.
# Complete the specification to do variance targeting
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(
model = "sGARCH",
variance.targeting = TRUE
),
distribution.model = "std"
)
# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
# Print the GARCH model implied long run volatility
sqrt(uncvariance(garchfit))## [1] 0.01099808# Verify that it equals the standard deviation (after rounding)
all.equal(sqrt(uncvariance(garchfit)), sd(sp500ret), tol = 1e-4)## [1] "Mean relative difference: 0.0003076497"# Specify model with AR(1) dynamics, GJR GARCH and skewed student t
flexgarchspec <- ugarchspec(
mean.model = list(armaOrder = c(1, 0)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "sstd"
)
# Estimate the model
flexgarchfit <- ugarchfit(data = sp500ret, spec = flexgarchspec)
# Complete and study the statistical significance of the estimated parameters
round(flexgarchfit@fit$matcoef, 6)## Estimate Std. Error t value Pr(>|t|)
## mu 0.000340 0.000083 4.114964 0.000039
## ar1 -0.029168 0.011651 -2.503390 0.012301
## omega 0.000001 0.000000 2.433509 0.014953
## alpha1 0.001521 0.003511 0.433308 0.664791
## beta1 0.918116 0.007553 121.549276 0.000000
## gamma1 0.143252 0.003340 42.884298 0.000000
## skew 0.917839 0.014330 64.048800 0.000000
## shape 6.932397 0.499364 13.882460 0.000000# str(flexgarchfit)A comment for the skew parameter: The skewed student t skewness parameter is statistically significantly different from zero. However 0 is not an interesting value to compare with for the skewness parameter. We should compare it with 1, since the the distribution is symmetric if the skewed student t skewness parameter equals one. Here the estimated skewness parameter is 0.917 and its standard error equals 0.0143. The estimate is thus more than two standard errors away from one and thus we can conclude that the skewness parameter is different from one and that the distribution is thus not symmetric.
# Specify model with constant mean, standard GARCH and student t
tgarchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "gjrGARCH", variance.targeting = TRUE),
distribution.model = "sstd"
)
# Fix the mu parameter at zero
setfixed(tgarchspec) <- list("mu" = 0)
# Estimate the model
tgarchfit <- ugarchfit(data = sp500ret, spec = tgarchspec)
# Verify that the differences in volatility are small
plot(sigma(tgarchfit) - sigma(flexgarchfit))
Comparison of mean errors and variance errors:
# Compute prediction errors
garcherrors <- residuals(garchfit)
gjrerrors <- residuals(tgarchfit)
# Compute MSE for variance prediction of garchfit model
mean((sigma(garchfit)^2 - garcherrors^2)^2)## [1] 1.293864e-07# Compute MSE for variance prediction of gjrfit model
mean((sigma(tgarchfit)^2 - gjrerrors^2)^2)## [1] 1.216576e-07# Print the number of estimated parameters
print(coef(garchfit))## mu alpha1 beta1 shape omega
## 6.732457e-04 7.108829e-02 9.232537e-01 6.377146e+00 6.843804e-07print(coef(tgarchfit))## mu alpha1 beta1 gamma1 skew shape
## 0.000000e+00 1.556731e-06 9.154120e-01 1.533007e-01 9.029445e-01 7.218296e+00
## omega
## 1.311604e-06or:
# Print the number of estimated parameters
length(coef(garchfit))## [1] 5length(coef(tgarchfit))## [1] 7# Print likelihood of the two models
likelihood(garchfit)## [1] 24271.46likelihood(tgarchfit)## [1] 24366.18# Print the information criteria of the two models
infocriteria(garchfit)##
## Akaike -6.642251
## Bayes -6.638475
## Shibata -6.642251
## Hannan-Quinn -6.640952infocriteria(tgarchfit)##
## Akaike -6.667901
## Bayes -6.663182
## Shibata -6.667902
## Hannan-Quinn -6.666279Diagnose the standardize absolute returns
# Compute the standardized residuals
stdret <- residuals(tgarchfit, standardize = TRUE)
# Compute their sample mean and standard deviation
mean(stdret)## [1] 0.03628148sd(stdret)## [1] 1.004238# Compute the standardized returns
stdsp500ret <- residuals(tgarchfit, standardize = TRUE)
# Compute their sample mean and standard deviation
mean(stdsp500ret)## [1] 0.03628148sd(stdsp500ret)## [1] 1.004238# Correlogram of the absolute (standardized) returns
par(mfrow = c(1, 2))
acf(abs(sp500ret), 22)
acf(abs(stdsp500ret), 22)
# Ljung-Box test
Box.test(abs(stdsp500ret), 22, type = "Ljung-Box")##
## Box-Ljung test
##
## data: abs(stdsp500ret)
## X-squared = 39.116, df = 22, p-value = 0.01369Shorter sample:
# Estimate the model on the last 1000 observations
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "sGARCH"),
distribution.model = "std"
)
garchfit1000 <- ugarchfit(data = tail(sp500ret, 1000), spec = garchspec)
round(garchfit1000@fit$matcoef, 4)## Estimate Std. Error t value Pr(>|t|)
## mu 0.0007 0.0002 4.2618 0.000
## omega 0.0000 0.0000 1.2107 0.226
## alpha1 0.2138 0.0478 4.4739 0.000
## beta1 0.7670 0.0453 16.9136 0.000
## shape 4.7618 0.7632 6.2395 0.000# Compute standardized returns
stdret <- residuals(garchfit1000, standardize = TRUE)
# Do the Ljung-Box test on the absolute standardized returns
Box.test(abs(stdret), 22, type = "Ljung-Box")##
## Box-Ljung test
##
## data: abs(stdret)
## X-squared = 22.99, df = 22, p-value = 0.4023In three steps:
# Estimate the GARCH model using all the returns and compute the in-sample estimates of volatility
garchinsample <- ugarchfit(data = sp500ret, spec = garchspec)
garchvolinsample <- sigma(garchinsample)
# Use ugarchroll for rolling estimation of the GARCH model
garchroll <- ugarchroll(garchspec,
data = sp500ret,
n.start = 2000, refit.window = "moving", refit.every = 2500
)
# Set preds to the data frame with rolling predictions
preds <- as.data.frame(garchroll)
# Compare in-sample and rolling sample volatility in one plot
garchvolroll <- xts(preds$Sigma, order.by = as.Date(rownames(preds)))
volplot <- plot(garchvolinsample, col = "darkgrey", ylim = c(-0.1, 0.1))
volplot <- addSeries(garchvolroll, on = 1, col = "blue")
plot(volplot, main = "In-sample versus rolling vol forecasts")
# ? not visibleThe proof of the pudding is in the eating.
# Inspect the first three rows of the dataframe with out of sample predictions
garchpreds <- as.data.frame(garchroll)
head(garchpreds, 3)## Mu Sigma Skew Shape Shape(GIG) Realized
## 1996-11-29 0.0006421517 0.006061525 0 5.238311 0 0.0026754967
## 1996-12-02 0.0006421517 0.006000423 0 5.238311 0 -0.0006076458
## 1996-12-03 0.0006421517 0.005935000 0 5.238311 0 -0.0109442741# Compute prediction errors
e <- garchpreds$Realized - garchpreds$Mu
d <- e^2 - garchpreds$Sigma^2
# Compute MSE for the garchroll variance prediction
garchMSE <- mean(d^2)
# Compute MSE for gjrgarchroll
gjrgarchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "std"
)
gjrgarchroll <- ugarchroll(gjrgarchspec,
data = sp500ret,
n.start = 2000, refit.window = "moving", refit.every = 2500
)
gjrgarchpreds <- as.data.frame(gjrgarchroll)
e <- gjrgarchpreds$Realized - gjrgarchpreds$Mu
d <- e^2 - gjrgarchpreds$Sigma^2
gjrgarchMSE <- mean(d^2)
(gjrgarchMSE - garchMSE) / garchMSE## [1] -0.04899196# When doing rolling sample GARCH model predictions, the comparison is fair. You only use the returns available at the time of prediction. As such you are not exposed to look-ahead bias nor to overfitting, since the data used for the prediction is different from the actual return used to evaluate this prediction. The best predicting model wins. Depending on the application, the winner can be a simple model or a complex model.VaR:
# Extract the dataframe with predictions from the rolling GARCH estimation
garchpreds <- as.data.frame(garchroll)
# Extract the 5% VaR
garchVaR <- quantile(garchroll, probs = 0.05)
# Extract the volatility from garchpreds
garchvol <- xts(garchpreds$Sigma, order.by = time(garchVaR))# Analyze the comovement in a time series plot
garchplot <- plot(garchvol, ylim = c(-0.1, 0.1), col = "red")
garchplot <- addSeries(garchVaR, on = 1, col = "blue")
garchplot <- addSeries(sp500ret, on = 1, col = "grey")
plot(garchplot, main = "Daily vol and 5% VaR")
# Take a default specification a with a normal and skewed student t distribution
normgarchspec <- ugarchspec(distribution.model = "norm")
sstdgarchspec <- ugarchspec(distribution.model = "sstd")
# Do rolling estimation
normgarchroll <- ugarchroll(normgarchspec,
data = sp500ret,
n.start = 2500, refit.window = "moving", refit.every = 2000
)
sstdgarchroll <- ugarchroll(sstdgarchspec,
data = sp500ret,
n.start = 2500, refit.window = "moving", refit.every = 2000
)
# Compute the 5% value at risk
normgarchVaR <- quantile(normgarchroll, probs = 0.05)
sstdgarchVaR <- quantile(sstdgarchroll, probs = 0.05)
# Compute the coverage
actual <- xts(as.data.frame(normgarchroll)$Realized, time(normgarchVaR))
mean(actual < normgarchVaR)## [1] 0.05658415mean(actual < sstdgarchVaR)## [1] 0.05887248# The coverage is clearly better when using a skewed student t distribution. Modelling correctly the distribution is of utmost important for getting the quantiles (and thus the value-at-risk) right. Note that in the implementation we took 2000 observations to reduce the computation time in this example. In practice you will be able to improve performance by re-estimating at a higher frequency.Fit the model to a particular period:
# Estimate the model
garchfit <- ugarchfit(data = sp500ret["/2006/12"], spec = garchspec)
# Fix the parameters
progarchspec <- garchspec
setfixed(progarchspec) <- as.list(coef(garchfit))
# Use ugarchfilter to plot the estimated volatility for the complete period
garchfilter <- ugarchfilter(data = sp500ret, spec = progarchspec)
plot(sigma(garchfilter))
Compare the forecasts for different periods
garchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(
model = "sGARCH",
variance.targeting = FALSE
),
distribution.model = "std"
)
# Compare the 252 days ahead forecasts made at the end of September 2008 and September 2017
progarchspec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0)),
variance.model = list(model = "sGARCH"),
distribution.model = "std"
)
### check how it works now:
garchfit2008 <- ugarchfit(data = sp500ret["/2008-09"], spec = garchspec)
garchfor2008 <- ugarchforecast(fitORspec = garchfit2008, n.ahead = 252)
garchfit2017 <- ugarchfit(data = sp500ret["/2017-09"], spec = garchspec)
garchfor2017 <- ugarchforecast(fitORspec = garchfit2017, n.ahead = 252)
# plot
par(mfrow = c(2, 1), mar = c(3, 2, 3, 2))
plot(sigma(garchfor2008), main = "/2008-09", type = "l")
plot(sigma(garchfor2017), main = "/2017-09", type = "l")
# we apply here an earlier specification with GARCH-in-Mean gim_
gim_garchfit <- ugarchfit(data = sp500ret["/2008-09"], spec = gim_garchspec)
garchforecast2008 <- ugarchforecast(data = sp500ret["/2008-09"], fitORspec = gim_garchfit, n.ahead = 252)
gim_garchfit <- ugarchfit(data = sp500ret["/2017-09"], spec = gim_garchspec)
garchforecast2017 <- ugarchforecast(data = sp500ret["/2017-09"], fitORspec = gim_garchfit, n.ahead = 252)
# fitORspec = garchspec
par(mfrow = c(2, 1), mar = c(3, 2, 3, 2))
plot(sigma(garchforecast2008), main = "/2008-09", type = "l")
plot(sigma(garchforecast2017), main = "/2017-09", type = "l")
Simulations
# A code to simulate 5 time series of 10 years of daily returns
spec <- ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "sstd",
fixed.pars = list(mu = 0.001, omega = 0.00001, alpha1 = 0.05, beta1 = 0.90, shape = 8, skew = 2)
)
# simulate the path
# path.sgarch = ugarchpath(spec, n.sim=3000, n.start=1, m.sim=1)
simgarch <- ugarchpath(spec = spec, m.sim = 5, n.sim = 10 * 252, rseed = 210)
# Plot the simulated returns and volatility of the four series
simret <- fitted(simgarch)
plot.zoo(simret)
plot.zoo(sigma(simgarch))
# Compute the corresponding simulated prices and plot them
simprices <- exp(apply(simret, 2, "cumsum"))
matplot(simprices, type = "l", lwd = 3)
# Estimate using default starting values
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
# Print the estimated parameters and the likelihood
coef(garchfit)## mu omega alpha1 beta1 shape
## 6.751676e-04 6.167870e-07 7.470453e-02 9.230754e-01 6.087123e+00likelihood(garchfit)## [1] 24272.28# Set other starting values and re-estimate
setstart(garchspec) <- list(alpha1 = 0.05, beta1 = 0.9, shape = 6)
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)
# Print the estimated parameters and the likelihood
coef(garchfit)## mu omega alpha1 beta1 shape
## 6.706921e-04 6.143096e-07 7.436183e-02 9.234476e-01 6.110132e+00likelihood(garchfit)## [1] 24272.27#Covariance - with the dataset
# Compute the standardized returns, together with their correlation
stdret_aparch <- residuals(aparchfit, standardize = TRUE)
gim_garchfit <- ugarchfit(data = sp500ret, spec = gim_garchspec)
stdret_gim <- residuals(gim_garchfit, standardize = TRUE)
garchfit_gim_cor <- as.numeric(cor(stdret_aparch, stdret_gim))
print(garchfit_gim_cor)## [1] 0.9912919Other series - home:
-download data for MSFT and WMT -fit models -show correlations in the moving window